Hydrodynamic behaviour of an Abelian Sandpile Model with Laplacian rules 
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We present a sandpile model, in which the instability of a site is determined also by the variables 
in a neighbourhood. This is a modification of the Abelian Sandpile Model, in which abelianity is 
preserved: it shares several mathematical properties of the original abelian model, while producing 
a more realistic dynamics. We show how our model presents interesting hydrodynamic features. 



I. INTRODUCTION 

Avalanche processes describe a variety of phenomena in 
Nature, ranging from, of course, avalanches (e.g. in piles 
of sand), to epidemic spread, allometric growth, . . . They 
are considered as a possible theoretical basis of 1/ / noise 
(pink noise), and of the appearence of Self-Organised 
Criticality in physical and natural systems [1, 2]. 

A lattice modelisation of non-equilibrium avalanche 
dynamics is expected to be based on local toppling and 
diffusion rules: if a toppling threshold inequality is sat- 
isfied at x (i.e., x is unstable, e.g., the slope at x is too 
steep), a diffusion process occurs at x, which may trigger 
further topplings, and so on, producing an avalanche. 

When more sites are simultaneously unstable, an or- 
dering prescription must be given for definiteness. This 
produces a variety of (stochastic and deterministic) pos- 
sible rules, that share the same crucial classification dif- 
ficulties of non-equilibrium Statistical Mechanics, in con- 
trast with the more clear and rigorous understanding of 
universality in equilibrium Critical Phenomena. 

A simple dynamics, introduced by Bak, Tang and 
Wiesenfeld [1], involves a fixed, and ultra-local thresh- 
old rule, of the form z(x) > h (instead of a Laplacian 
one, that involves heights z{y) at neighbours y ~ x). Re- 
markably, as shown in the work of Dhar and collaborators 
[3,4], for such a dynamics topplings at different sites com- 
mute, this justifying the name of Abelian Sandpile Model 
(ASM) given to this class of systems. The ingredient of 
abelianity is crucial in the determination of several hid- 
den mathematical features of the model, among which, 
notably, a bijection with spanning trees, configurations in 
an equilibrium Statistical Mechanics model. Abelianity 
also allows to achieve scaling to the continuum limit, ei- 
ther by the forementioned correspondence, or by the fact 
that adding particles and relaxing in rounds is equivalent 
to add all the particles at once, and perform a unique 
complete relaxation, this latter process being easier to 
analyse theoretically in the large- volume limit (besides 
that algorithmically convenient). 

Unfortunately, abelianity seems to be a fragile ingre- 
dient. Small modifications of threshold or diffusion rules 
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easily break commutativity. In particular, except for a 
construction specific to 1-dimcnsional chains, it seems 
impossible to implement genuine and realistic toppling 
rules, based on gradient or Laplacian, within ASM's. In 
[5] we presented a new class of ASM's, allowing for multi- 
site topplings: this increases the spectrum of possibilities, 
but still within ultra-local threshold rules. 

Here we present a family of ASM's based on Laplacian- 
like toppling rules, possibly non-linear. The most natu- 
ral rule should be linear homogeneous, e.g. of the form 
z(x) > {z(y)) yr ^ x . However we need to introduce non- 
linearity, both for preserving abelianity, and for having a 
compact space of configurations. The simplest variant is 
of the form z(x) > a (z{y)) yr ^ x + S, thus linear inhomo- 
geneous. Other variants are discussed later on. 

Furthermore, in order to preserve abelianity, we need 
to adopt non-compact diffusion rules: if ordinary nearest- 
neighbour diffusion is (Df)(x) = J2 y ~ x (f(v) ~ f( x ))> 
we need a more general diffusion operator, of the form 
(D u f)(x) = J2 y u(x - y)(f(y) - f(x)). Fortunately, 
finite-range functions u(r), e.g. u(r) ~ exp(— A|r|), are 
allowed, and should show the same phenomenology of 
nearest-neighbour ones. 

Interestingly, even in our largest class of models, con- 
figurations recurrent under a dynamics of random in- 
creases of the heights (plus relaxation), as in ordinary 
ASM's, still have an uniform steady-state probability dis- 
tribution, and a generating function given by the Kirch- 
hoff matrix-tree formula. It is possible that natural bijec- 
tions exist with spanning trees also in our wider context. 

Thus we have continuous- variable sandpiles with fixed- 
threshold (F) or non-local Laplacian-like (L) rules, and 
with compact (C) or non-compact (N) diffusion. Of the 
four possibilities, the three (FC), (FN) and (LN) have 
abelian realisations, with (FC) being the "ordinary" case. 
We will use the shortcuts X-, Y-, XY-ASM (X=F,L; 
Y=C,N) for these classes of models, and ASM for the 
original model with discrete variables. 

In order to support our claim that Laplacian-likc 
threshold rules are the crucial new ingredient here (while 
continuous variables and non-compact diffusion are a 
technical accident), we show how LN-ASM models have 
interesting phenomcnological properties when observed 
under "hydrodynamic" experiments, that are not shared 
neither with the ordinary ASM and FC-ASM, nor with 
the FN-ASM models. 
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II. THE MODEL 

For simplicity we will only analyse translationally invari- 
ant sandpiles with continuum height variables, on a por- 
tion V = < S>i< a <d'^' L c °f tne ^-dimensional hypercubic 
lattice, with periodic or open boundary conditions. Ex- 
tensions, to be discussed elsewhere, could include models 
defined on arbitrary graphs, even directed, and realisa- 
tions with discrete variables. 

We start by considering dynamics in which the heights 
are increased, stochastically or deterministically, and 
avalanches are possibly produced. When we say that a 
set of configurations is "left stable by the dynamics" , we 
mean by any dynamics of this form. In particular, this 
excludes the inverse-toppling dynamics discussed in [5] . 

We denote by w(£) the Laplace transform of w(x). 
We say that a function f(x) : Z d — > E is symmetric if it 
has the symmetries of the cubic lattice, f(x\, . . . , x d ) = 
f(e 1 x a(1) ,. . .,e d x a(d) ), for a E & d and e G {±l} d 

Our sandpile model is determined by two non-negative 
symmetric functions w(x) and u(x), with w(0) = u(0) = 
0, a dissipation parameter /i > 0, and two functions f(z), 
g(z), strictly- and weakly-monotonic respectively, with 
g(z) having a finite Lipschitz constant £. We say that x 
is unstable if 



f(z(x)) > w (v) 9{z{x + y)) . 



(1) 



In such a case, a toppling t x may occur at x, modifying 
the height function as 

t x • z(y) -»• | z{x) _ {1 + fl) J2 yl u (y>) y = x ^ 

Note that \x > implies that J2 X z(x) strictly decreases 
at each toppling. A configuration z(x) is stable if no x is 
unstable. We have the F-ASM if f(z) = z and g(z) = 1. 
Furthermore, it is customary to take w(x) = u{x) = 
8\ x \_i the indicator functions on nearest neighbours, this 
producing a FC-ASM. We have an obvious affine covari- 
ance z — > jiz + 70, with 71 > 0, and invariance under 
multiplication of the threshold rule (1) by a positive con- 
stant, that we exploit later on. 

We require the sandpile to have three properties: 

[A] A positive cone ft — {z : z(x) > z m \ D for all x} is 

left stable by the dynamics. 

[B] The set of stable configurations within fl has finite 

non-zero measure. 

[C] Within fi, the topplings are abelian. 

Call C = (1 + M ) E y ^ <V) > and w(0) = £ y w(y). 
We can use the covariance to fix z m ; n = 0, i.e. O = . 
In this case, condition [A] means that, if z(y) > for all 
y, and x is unstable, then (1) implies z(x) > C, that is, 
using the monotonicity, C < C := f' 1 (g(0)w(0)) . 



One easily sees that the set of stable configurations 
has non-zero measure, because (0, C'] v is a subset. Sup- 
pose that the limit for h — > +00 of f(h) — g(h)w(0) 
exists. If it is positive, there exists h max such that 
f(h) > g(h)w(0) for all h > ft- max , and any configura- 
tion in n with max y z(y) > ft. max is unstable at least 
at the position of the max, so stable configurations are 
contained within (0,ft. max ] y and [B] is verified. If it is 
negative, there exists e > 0, and an unbounded set of h's 
in M+, such that the cubes h — e < z(y) < h are sta- 
ble, because f(h) < (g(h) — e£)w(0) < g{h — e)w(O), and 
[B] is not verified. When the limit is equal to zero, or 
does not exist, the condition needs to be analysed more 
deeply, and we don't do this here. 

A sufficient condition for [C], using the result on the 
analysis of [A] and the Lipschitzianity of g(z), is that, 
for all x and all z > C , 

f(z + u(x))-f(z)>£^2w(y)u(x-y). (3) 

Assume that f'(z) > £' > for all z > C . Then we have 
the condition, for all x ^ 0, 



£' u(x) > £{w *u) (x) 



(4) 



It is easy to see that, if u{x) = uq cxp(— A | as | ) , then 
sup^^p (w*u){x)/u(x) — w(Xn) for n some unit vector. 
For such a function u the condition reads w{\n) < £' /£. 
From now on, we will only consider u's of this form. 
Clearly w(Xn) > w(0) for all A > 0, and the difference 
goes to zero in the limit, thus some value A > satisfying 
the condition exists if and only if w(0) < £' /£. 

Use affine covariance to set C = 1. We can eliminate 
w(0) to get the sufficient and necessary condition on / 
and g, to admit w, A producing a LN-ASM 
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lim 



(5) 



Recall that £' = min z> c f'(z) and £ = max z>0 g'{z). 
One easily sees /(M) > Ml' + O(l) and g(M) < 
M£ + 0(1), implying that the minimum above is always 
realised on the first quantity, i.e. the condition [C] is al- 
ways stronger than [B] . In the case of a linear threshold 
rule, g(z) = g + zgi and f(z) = f + zfx (with g , fi > 
and <7i > 0) we get the constraint (/o + /i)<?i < /i<?o- 
If f(z) is quadratic, with / 2 > 0, we get analogously 

(/o + /i+/ 2 )3i < (A + 2/2)50. 

We say that a realisation of sandpile model as above is 
tight if conditions [A] and [C] are satisfied in a tight way. 
Assuming that f'(z) is monotonic, that g(z) = g + zg\ 
is linear, and setting C = 1, this means that w and A are 
set to satisfy 
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= w(0) < w(Xn) 
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(6) 



For W C V, define xw{x) = 1 if x e W and oth- 
erwise. It is easy to see that the frame identity Id/ 
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of the LN-ASM on V associated to (/, g, w, u) depends 
only on u and the domain V (in particular, it is the 
same as in the associated FN- ASM), and is given by 
Id/(x) = — (xv * u ){ x ), the difference in height after 
one toppling has occurred at each site. Conversely, the 
recurrent identity Id r depends also on w, / and g. 

III. MULTIPLE THRESHOLD RULES 
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For fixed w and u, we may have families of LN-ASM, 
which provide a continuous deformation of the FN- ASM 
with threshold rule z{x) > h, and thus f(z) — z and 
g(z) = h/w(0). 

We may extend the formalism of the previous section to 
functions /, g : R — > R s , for s > 1, and, in (1), replace > 
with y, the canonical partial ordering (/ y g iff f a > g a 
for all 1 < a < s). This provides a different mechanism 
for continuous deformations of F-ASM's: we can take s > 
2, and keep the original rule for one of the components. 

In the analysis of properties [A] and [C], the dif- 
ferent components have a non-trivial interplay in a 
unique respect. For property [A] we need C < C' a :— 
/~ 1 (g Q (0)w(O)), for all a. Then, for property [C] we 
need w(\n) < £' a /£ a , for all a. Apparently, these con- 
straints are factorised, but this is not the case. While £ a , 
the Lipschitz constant of g a {z), is in fact independent 
from the other components, l' a is defined in terms of the 
allowed range for unstable heights, more precisely 



*«■■= , J* (£(*))■ 

0>max,3(G^) 



(7) 



Say that the (one or more) indices realising the max of 
Ci are the leading components. For a valid realisation 
of the sandpilc in dimension s > 2, the restriction of 
the threshold equation to a subset of components con- 
taining at least one leading component still produces a 
valid realisation. If we require non-redundancy, we get 
the constraint 

[D] For each 1 < a < s, there exists z £ SI such that 



fp(z(x)) - T, v w(y)9p(z(y)) < only for f3 



a. 



This constraint is specially important if we have s = 2, 
a = 1 is the fixed-threshold constraint z(x) > h, and is 
also the only leading component. Our sandpile model is 
genuinely distinct from the original version if and only 
if the condition is verified, and the point at which this 
happens in a tight way is the starting point of the family 
of continuous deformations. 



IV. ANALYSIS OF EXAMPLES IN d =1 AND 2 

The tight realisations of the sandpile on an infinite linear 
chain, with w(x) = u><5|.,,| ! and u(x) = uA~\ x \, require 
A > 1 and give 



A 2 + 1 _ f(l)g 
/(IK ' 



A-l 



w 



2ffo 



(8) 



FIG. 1: Time evolution of a centrally-seeded sandpile in d = 1, 
with random low-height initialisation. Right: in a LN-ASM 
as described in Setion IV, with A = 2. Left: in the associated 
FN- ASM. In all our figures, black/white = low/high. 



If also f(z) is linear, then we can set 

f(*) = *-^f?T-> 9(z)=z + l, (9) 



A 2 + l ' 



2A 



A-l 



w 



A 2 + l ' 



(10) 



The maximal possible height is h max = (A 2 + 1)/(A — 1 ) 2 . 

Figure 1 compares the time trace of a simple protocol 
(steady injection of sand in the middle of the lattice), 
between the model above, for A = 2, and the FN-ASM 
using the same function u(x). In both cases we observe 
a linear growth of a plateau. While in the F-ASM the 
heights in the plateau show fuzzy short-range fluctua- 
tions, in the L-ASM the profile is smooth, with no need 
of coarse-graining. 

In the two-dimensional case, we can still choose the 
model to be tight, w(x) = io5i x m, g(z) = z + 1, 
f(z) = z + /o, u(x) = wcxp(— X\x\), and C = C = 1 
(thus /o < 0, w = (l+/ )/4 and cosh(A) = £-1). As the 
series N(\) := X^x/o cxp(— A|a;|) has no closed form, the 
parameter u has to be set numerically to u = A/"(A) _1 . 
For example, choosing /o = —0.15 gives u = 0.1164... 
This is the realisation presented in the example of Fig- 
ure 2, and later on in the following section. Again, we 
observe smoothness properties only in the LN-ASM case. 
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FIG. 2: Relaxation of z(x) = 32/(1 + \x\ 2 ) in the FN-ASM 
(left), and of z(x) = 196/(1 + |cc| 2 ) in the LN-ASM (right). 
On top: the profile of the middle row. 
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FIG. 3: Left: a portion of a typical outcome of toppling- 
antitoppling dynamics in the ASM. Right: the analogous out- 
come in our LN-ASM. 

V. EXTREME REGIMES 

We have seen above how the sandpiles with Laplacian-likc 
toppling rules, when investigated in regimes imitating re- 
alistic experimental settings ( "hydrodynamic regimes"), 
show smoothness properties not shared with the ordinary 
versions. These settings are characterised by the presence 
of empty regions, collecting the outcome of avalanches, 
and making them not percolating. In this section we 
discuss how the system behaves when it is driven to- 
wards extreme regimes, in which the avalanche size is 
regularised only by the finitcness of the volume. In this 
case, the new model preserves some resemblance with the 
intriguing combinatorial features of the ordinary ASM. 

A natural candidate is the study of the recurrent iden- 
tity, e.g. on a L x L portion of the square lattice, with 
open boundary condition. This configuration has a cen- 
tral role in the mathematical analysis of the model (it is 
the identity of the abelian group associated to recurrent 
configurations) [6]. It is obtained through a non-local 
procedure, as the fixed point of iterated addition of frame 
identities, and is very sensitive to the shape of the do- 
main and boundary conditions. In the ordinary discrete 
ASM, it shows convergence to self-similar limit shapes, 
involving fractal Sierpinskij-like structures. An appro- 
priate discussion is too complex to fit here, and we start 
instead with analysing its "building blocks", the strings 
[7]. These are one-dimensional structures in ASM's on 
bidimensional lattices, that also play a role in the clas- 
sification of the emerging periodic patterns under deter- 
ministic protocols [8, 9]. 

In the ordinary (discrete) ASM, strings first appear 
when performing "toppling-antitoppling" operations [5, 



7] on the maximally-filled configuration (z( max )(a;) = 3 
for all x, on the square lattice). For LN- ASM's we first 
need to generalise the concept of z( max )(:c), which turns 
out to be non-homogeneous. It is characterised as the 
unique solution to the system in which the threshold in- 
equalities (1) are replaced by equalities. Thus, it is sta- 
ble, but z (max) + ex{ x } is unstable for all x E V, e > 0. 

If /, g are linear and w is nearest-neighbour, the corre- 
sponding system has the form (A + a)z^ max '(a;) = b(x), 
plus vanishing boundary conditions. In our geometry, 
this is easily solved in Fourier basis and by method of 
images: the solution is smooth, and near to ft. max , the 
maximal possible height, for sites far from the boundary. 

Next, we need to define toppling-antitoppling opera- 
tors v x in continuum ASM's. Call z' the relaxation of 
z + sx{ x }, where s is an amount making x barely un- 
stable. Then call z" the antirelaxation of z' — s\{x}- 
We set v x z := z". The action of continuous toppling- 
antitoppling operators v x on z ( max ) ; i n LN-ASM, is anal- 
ogous to the action of discrete toppling-antitoppling op- 
erators a x a x in ASM. Examples are shown in Figure 3. 

The reason is that strings correspond to discontinuities 
in piecewise-linear toppling matrices, and even in contin- 
uous sandpiles the toppling matrices are discrete. This 
fact combines with the observations in [7, 9]. 

An aspect of this resemblance of ASM and LN-ASM 
occurs also in a conceptually-simpler protocol: add sand 
to an uniform configuration z{x) = ft. max , at a unique 
site, then relax. Typical outcomes are shown in Figure 4. 
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FIG. 4: Relaxing kL particles added at x — (L/3, L/4) to 
2(03) = ftmax, on a L x L square (here L — 81). Left: ASM, 
and k = 2. Right: GN-ASM, and k = 4. Graytones are scaled 
by z [m ^\x). 
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